Development and Numerical Analysis of 
"Black-box" Counterpropagating Wave Algorithm for 
Exact Quantum Scattering Calculations 
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In a recent series of papers [J. Chem. Phys. 121 4501 (2004), J. Chem. Phys. 124 034115 
(2006), J. Chem. Phys. 124 034116 (2006)] a bipolar counter-propagating wave decomposition, 
9 — 9+ + 9-, was presented for stationary bound states <& of the one-dimensional Schrodinger 
equation, such that the components 9± approach their semiclassical WKB analogs in the large 
action limit. The corresponding bipolar quantum trajectories are classical-like and well-behaved, 
even when 9 has many nodes, or is wildly oscillatory. In this paper, the earlier results are used to 
construct a universal "black-box" algorithm, numerically robust, stable and efficient, for computing 
accurate scattering quantities of any quantum dynamical system in one degree of freedom. 



I. INTRODUCTION 

In a recent series of three articles j^ 2 -! 3 - the use of bipo- 
lar counter-propagating wave methods (CPWMs)^^^^ 
for exact numerical solution of the time-independent 
Schrodinger equation was explored. The basic underly- 
ing idea is to decompose the stationary wavefunction, "J, 
into a sum of two counter- propagating wave components, 



(1) 



If done appropriately, such a decomposition can lead 
to very important ramifications for quantum trajectory 
methods (QTMs^i&^^^i^i.e., trajectory- 
based numerical techniques for performing exact 
quantum dynamics calculations, based on Bohmian 
mechanic o 15 i 16 ' 17 i 18 ' 19 i 20 — due to nonlinearity of the 
Bohmian equations of motion. In particular, the earlier 
series of articles has culminated in a set of trajectory- 
based time-dependent methods for computing station- 
ary scattering quantities (the theoretical underpinning of 
all chemical reactions) in one degree of freedom (DOF). 
These methods were found to be numerically stable and 
efficient for model systems within a certain regime of 
system parameters and accuracy typical of molecular 
applications. 2,3 However, more recent investigations have 
uncovered certain limitations that arise in more extreme 
cases. This has motivated the goal of the present pa- 
per, to perform a much more detailed numerical analysis 
of the methods, and to develop a single stable, robust, 
and efficient "black-box" algorithm that can be applied 
to virtually any 1 DOF system with a minimum of user 
intervention. 

In recent years, QTM's have arisen as a very 
promising tool for performing accurate quantum dy- 
namics calculations for many-DOF systems, using tra- 
jectory ensembles in a manner similar to classical 
simulations.™ In the conventional unipolar formula- 
tion of Bohmian mechanicS ) 15 i 16 ' 17 i 18 ' 19 i 20 the Madelung- 
Bohm amplitude-phase decomposition of the 1 DOF 



i.e. 



9(x,t) = R(x,t)e lS{x > t),h , 



(2) 



is used to generate the (presumably classical-like) quan- 
tum field functions, R(x,t) and S(x,t). For slowly- 
varying potentials, the classical field functions are also 
slowly- varying, except in the vicinity of caustics. The 
corresponding quantum field functions of Eq. ([2]) behave 
similarly when there is no interference. However, interfer- 
ence introduces non-classical-like oscillations in the quan- 
tum field functions, which in turn lead to numerical diffi- 
culties for QTM calculations — collectively referred to as 
"the node problem"^iS (even though true nodes per se 
need not be present). As interference is formally unavoid- 
able in any reactive scattering context, the node problem 
is to date the most important issue impeding the progress 
of QTM's as a general and robust tool for scattering dy- 
namics applications, although much progress has been 
made^^^ 

The bipolar CPWM approach to the node problem, 
as adopted here, and in the earlier series of articles ^i 3 . 
is inspired by classical and semiclassical theories — which 
employ a multipolar, Eq. ([lj-like decomposition as the 
means of preserving smooth and slowly-varying field 
functions. For stationary 1 DOF applications, the , J+ 
and "J- components are respectively, left- and right- 
traveling waves, with equal and opposite (group) veloci- 
ties, v(x) and —v{x). Oscillatory interference behavior is 
not evident in the individual component field functions, 
but arises naturally from their linear superposition. For 
stationary bound state calculations, cancellation of flux 
requires that the component densities, p+{x) = |\l/ + (a;)| 2 
and p-(x) = \^-(x)\ 2 , be equal. A reasonable gener- 
alization of classical field properties can be used to ar- 
rive at a unique exact quantum decomposition of the 
Eq. JT]) form, which — unlike the conventional unipolar 
ansatz [Eq. @] — satisfies the correspondence principle 
in the large-action limiti 

For stationary scattering, or continuum state calcula- 
tions, the non-uniqueness and non-square-integrability of 
the quantum solutions complicate matters somewhat,- 2,3 
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as might be expected. Boundary conditions play an es- 
sential role, and specifically for the usual left-incident 
boundary conditions presumed here, necessarily imply 
asymmetry and non-zero flux for the corresponding ex- 
act quantum stationary solution. Conventional scatter- 
ing theory already provides an asymptotic (i.e. |x| — > 
oo ) bipolar decomposition into the familiar "incident," 
"transmitted," and "reflected" plane wave components — 
in fact the only decomposition possible that does not 
result in asymptotically oscillatory field functions. An 
exact quantum bipolar CPWM decomposition that re- 
spects all of the above can be constructed^^ with 
joining the incident and transmitted waves continuously 
through the interaction region, and ^-(x), the reflected 
wave, damping to zero as x — ► oo. This is again achieved 
through analogy with semiclassical mechanics, albeit a 
"sophisticated" version . 25 i 26 i 27 A key feature is that the 
trajectories are actually classical, with all quantum ef- 
fects arising not from a quantum potential (as in conven- 
tional Bohmian mechanics), but through dynamical 1 §>± 
coupling. Conceptually, this gives rise to instantaneous 
trajectory hopping 28 from one CPWM component to the 
other, in the interaction region where the potential V(x) 
is changing. In effect, one has a local, time- dependent 
theory of stationary scattering, instead of the more tra- 
ditional global, time-independent picture. 

The new theory thus provides some pedagogical in- 
sight, analogous to ray optics and cavity ring-down 
experiments! 2 -^ However, this paper focuses primarily 
on numerical aspects of the new approach — which may 
be regarded as a relaxation method with exponentially- 
fast convergence, for computing stationary scattering 
states of desired boundary conditions, and associated re- 
active scattering quantities such as reflection and trans- 
mission probabilities. In comparison with other quan- 
tum scattering methods, the bipolar CPWM approach 
offers some decided numerical advantages. In particu- 
lar, there is no need to invoke complex scaling (analyt- 
ical continuation) 30 i 31 ' 32 or absorbing boundary condi- 
tions (optical potentials) f 3 ^^^^^^^ as trajectories 
and field quantities are real- valued throughout. More- 
over, the scaling of computational (CPU) effort is linear 
with the grid size, N, rather than proportional to N 3 . 

The new approach may therefore lead to some of the 
most efficient exact quantum scattering algorithms in ex- 
istence, but certain practical limitations must first be 
overcome. To begin with, the previous papers propose 
several distinct algorithms; a detailed numerical analysis 
must be performed to determine which of these is "best," 
from the perspective of numerical stability, robustness, 
and efficiency. In particular, some of these algorithms 
appear to be ineffective beyond a certain desired level of 
accuracy (Sec. Ill Bl) . The reasons for this must be elu- 
cidated via systematic study, in order to develop more 
efficient implementations. An even greater motivation, 
though, is the fact that the above problem is greatly ex- 
aggerated in the multidimensional case, resulting in nu- 
merical instabilities that have thus far hindered efforts in 



this direction. Another limitation is that some of the pre- 
vious algorithms are not applicable to completely general 
systems (Sec. Ill A"| . This is addressed here by generaliz- 
ing the methodology to allow for a nearly arbitrary choice 
of "classical" trajectory fSec. IIII Dl) . Finally, the robust- 
ness of the new algorithms must be demonstrated by ap- 
plication to representative molecular systems encompass- 
ing a very broad range of potentials, system parameters, 
and desired accuracy levels. One such application area, 
known to cause difficulties for conventional semiclassical 
and exact quantum methods, is particularly important — 
i.e., the deep tunneling regime (Sec. IIVB|) . Other exam- 
ples considered here include asymmetric potentials with 
barriers (Sec. IIVD|) and systems with reaction interme- 
diates fSec HVEl) . 



II. BACKGROUND 

A. Exact Quantum Dynamics Using 
Counterpropagating Trajectories 

Here, we summarize some of the developments of 
Refs. 1-3. The basic goal is to define a set of time- 
evolution equations for the two 1 DOF counterpropagat- 
ing waves, ty±(x,t), such that in the large t limit, the 
total ^>(x,t) of Eq. {1} approaches the exact stationary 
state evolution for the desired energy, E, and boundary 
conditions. This general approach is thus a "relaxation 
method," for which ty(x,t) at intermediate times need 
not adhere to the actual time-dependent Schrodinger 
equation — although there are versions of the method for 
which the latter property is also satisfied in a sense^i 3 - For 
the usual left-incident boundary condition, the initial re- 
flected wave must be zero, i.e. fy-(x,0) = 0. The initial 
ty+(x, 0) is essentially arbitrary, as it is presumed 2 -! 3 - that 
any reasonable choice will converge exponentially quickly 
in the large time limit. However, two natural choices 

have been considered: (1) \I> + (x, 0) 



exp 



i\/2mEx/h 

i.e. incident asymptotic plane wave extended through- 
out space [with V(x) — > presumed as x — > — oo]; 

(2) ^ + (x,0) = Q(xl — x) exp i\/2mEx/h~ , i.e. inci- 
dent wave truncated at left edge of interaction region, 
xl < x < xr. In Ref. H, (1) and (2) above are re- 
ferred to respectively as "non-wave-front" and "wave- 
front" versions. The latter gives rise to a useful cav- 
ity ring-down dynamical interpretation* 2 -^ although the 
former is found to converge more quickly to the exact so- 
lution, and is therefore used throughout this paper. Note 
that (1) also has its own interpretation, to be discussed 
shortly. 

From Ref. 0, we also learn that there are two dis- 
tinct semiclassical formalisms that can be applied to 
generate exact Eq. (Q} decompositions (and associ- 
ated time-evolution equations) with the desired theoret- 
ical properties — the Bremmer approach (B) 2 ^ and the 
Froman approach (F)i2S Of the two, the B approach can 
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only be applied when the trajectories are classical trajec- 
tories. Unfortunately, classical trajectories cause numer- 
ical instabilities for both B and F evolutions whenever 
there is tunneling, i.e. for all below-barrier energies, E, 
and is therefore not a viable choice for a completely "ro- 
bust" algorithm. The F approach, on the other hand, can 
be generalized for almost any desired choice of "semiclas- 
sical" trajectory. In particular, constant velocity trajec- 
tories (where the speed corresponds to the incident plane 
wave energy E) handle tunneling without difficulty, and 
lead to smoother field functions and remarkably simple 
F time-evolution equations: 

^ = i (£ - m± -iy* T (3) 

Some comments on Eq. (J3J) are in order. First, these 
are the Lagrangian, or "hydrodynamic" equations of mo- 
tion, meaning that the left hand-side is a total time 
derivative, as would be appropriate for a trajectory 
(moving-grid) calculation. Note that they require no dif- 
ferentiation of any kind — not even of the potential itself, 
to construct force fields. Note also that the final, cou- 
pling term, i.e. that responsible for all quantum correc- 
tions beyond the "semiclassical" approximation, is pro- 
portional to V(x), and vanishes in the x — > — oo reactant 
asymptotic limit. It also vanishes in the product asymp- 
tote provided the potential is "asymptotically symmet- 
ric," i.e. V(x) — >= as x — > oo (the asymmetric case 
will be addressed in Sec. IIII Cp . This ensures that the 
ty±(x) solutions exhibit the requisite plane wave behav- 
ior for the desired energy E. Finally, we note that for 
the special case V = 0, i.e. the potential that would 
classically give rise to the constant velocity trajectories 
used, all coupling vanishes. This implies that the non- 
wave-front initial condition (1) described above is in fact 
the "basic WKB"^2& approximation in the generalized F 
sense, as will be exploited in Sec. IIII Dl 

Since the ^ + and \&_ trajectories move in opposite 
directions, the two moving grids are necessarily incom- 
mensurate at an arbitrary time, which numerically ne- 
cessitates the use of interpolation for evaluating the cou- 
pling term in Eq. ([3]). Alternatively, an "Eulerian" or 
fixed-grid implementation can be adopted, for which the 
corresponding time-evolution equations are 

—± = Tv ^' ± + -(E-V)* ± - (4) 

where v — y^2E/m is the trajectory speed, and primes 
denote spatial differentiation. The grids are now com- 
mensurate at all times, obviating the need for interpo- 
lation, but numerical differentiation is now required to 
evaluate the new convective term in Eq. ((4]). 

B. Numerical Analysis of Fixed-grid Scheme 

Of the various possible numerical schemes presented in 
Ref. H and above, we have thus far settled on just two can- 



didates for the present purpose: the trajectory and fixed- 
grid versions of the non- wave- front, F, constant-velocity 
trajectory method. To choose one over the other, we 
must consider the stated criteria of numerical stability, 
efficiency, and robustness (i.e. generality of applicabil- 
ity across different systems and computational accuracy 
levels) . 

In Ref. [U, all of the calculations were applied to one 
degree-of- freedom (1 DOF) systems with similar param- 
eter values and energies, and for the most part, to a de- 
sired accuracy level of 10~ 4 . In this context, the two 
numerical schemes were found to be of roughly compa- 
rable efficiency (i.e. total CPU time), with the trajec- 
tory approach requiring substantially fewer grid points 
(around 4x per DOF) and larger time steps (x 10), but a 
much longer CPU time per step, due primarily to the 
nearest-neighbor search algorithm that was employed. 
On the other hand, there were also clear indications that 
the trajectory method would become superior at higher 
accuracies — e.g., Ref. H Fig. 5, which suggests that the 
fixed-grid error for reflection and transmission probabili- 
ties cannot be reduced below 10 or so, no matter how 
many grid points are used. 

A limit on achievable accuracy is clearly not a desirable 
trait for a robust method, and in particular, is complete 
anathema for accurate deep tunneling applications. Con- 
sequently, a detailed numerical analysis of this situation 
was performed for this paper, in order to gain a thorough 
understanding of the factors responsible. In the course 
of this, it was discovered that the limiting accuracy can 
be substantially worse for other reasonable system pa- 
rameter choices, especially for multidimensional systems, 
and/or those with reaction intermediates — lending ad- 
ditional imperative to the exercise. The manifestation 
of the error was quickly identified as oscillations in the 
right-asymptotic density p+(x) = |5' + (a;)| 2 , which form 
around x = xr, and then propagate inward over time — 
whereas it is clear from Eq. that the solution p' + (x) 
should be zero asymptotically. 

Ref. [H, and the above-mentioned preliminary numer- 
ical investigations for this paper, computed fixed-grid 
spatial derivatives using centered, fourth-order finite 
differences 4 ^ (one-sided for the endpoints xr and xl), 
and fourth-order Runge-Kutta time integration. 41 How- 
ever, the bulk of the present analysis uses second-order fi- 
nite difference (two-sided for interior grid points and one- 
sided for endpoints) and first-order forward Euler time 
integration, as a deliberate means of exacerbating the 
problem and simplifying the analysis, in order to identify 
the cause. The oscillation amplitudes were indeed found 
to be much larger in this case. It was also found that 
the oscillation wavelength was always 2-3 grid points, re- 
gardless of the density of grid points used. With respect 
to number of grid points, N, (keeping xl and xr fixed), 
the amplitude oscillation was found to decrease roughly 
inversely. Reducing the time step had no effect what- 
soever on the oscillation amplitude or wavelength, but 
did reduce the inward propagation speed, roughly pro- 
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portionally. 

All of the above findings support the notion that the 
oscillation problem is somehow due to spatial differenti- 
ation, ft might also be theorized that oscillations arise 
from the ^± phase difference inherent in the coupling 
term of Eq. (|3|), but the wavelengths are incommensu- 
rate with what is actually observed. To definitively rule 
out coupling as the source of the problem, however, a 
calculation was performed for the V = free particle sys- 
tem, for which there is no coupling — yet the above oscil- 
lations were still observed. Another alternative explana- 
tion would be von-Neumann instabilities. 42 A straightfor- 
ward theoretical analysis is possible for the free particle 
case, and so this was conducted, and compared to numer- 
ical tests for a wide range of time steps and grid spacings. 
In the theoretically von Neumann stable regime, the os- 
cillations were still observed, whereas in the von Neu- 
mann unstable regime, the numerical solutions diverged 
over time as expected — a different phenomenon entirely. 

In fact, none of the vast literature on PDE simulations 
of which the author is aware can account for the odd os- 
cillatory behavior observed. The only sensible explana- 
tion of the cause can therefore be the "mixed" boundary 
conditions, which are highly non-standard, and perhaps 
have not been considered previously. In particular, 



¥-(x Rt t)=0 



ty + (x L ,t) = exp 



-~n m 



(5) 



so that the two different components have boundary 
conditions on opposite sides. This situation is com- 
pletely natural, given the counterpropagating nature of 
the method, and cannot be changed through any straight- 
forward theoretical recasting. Within the fixed-grid 
framework, various numerical solutions have been es- 
sayed, including asymptotic plane wave rescaling, and 
amplitude/phase decomposition of the endpoint \&± com- 
ponents, but none appears to provide substantial relief 
from the oscillation problem described above. 



III. THEORY AND ALGORITHM 
DEVELOPMENT 

A. Constant-velocity Trajectory Scheme 

We are thus led naturally to the trajectory, rather than 
fixed-grid, version of the non- wave- front, F, constant- 
velocity method. Ref. Fig. 5 indicates the high ac- 
curacies (around 10~ 8 ) of which the trajectory approach 
is capable. This can be largely attributed to the fact that 
the trajectory version does not require numerical spatial 
differentiation. However, a more detailed numerical anal- 
ysis reveals that the errors in this case also manifest as 
asymptotic oscillations, albeit of much smaller amplitude 
than the fixed-grid case, and with different scaling. More 
specifically, for a test Eckart system comparable to that 
of Ref. H an inverse quartic relation is found between 
the oscillation amplitude and grid size, N. This can also 



be roughly inferred from Ref. Fig. 5, and is in any 
case not surprising, given that third-order interpolation 
is used (Sec. EH]). 

As in the fixed-grid case, the oscillation amplitudes and 
wavelengths are completely unaffected by time-step size. 
However, the latter quantity does turn out to be closely 
related to another, more serious type of numerical er- 
ror also observed: in the large t limit, the magnitude of 
^f^-(x,t) grows linearly with t in the left-asymptotic re- 
gion. For typical A values and desired accuracy levels, 
using fourth-order time evolution as in Ref. 0, the growth 
rate is so slow as to contribute insignificantly to the final 
error. However, this asymptotic growth does present a 
very significant challenge for extremely accurate calcu- 
lations (e.g. 10 -12 as in Sec. IIVB[) and must therefore 
be addressed. Once again, to exaggerate the effect for 
analytical purposes, first-order forward Euler time inte- 
gration was used. For the test Eckart system, the growth 
rate (with respect to system time, not number of time 
steps) was found to depend roughly linearly on A, and 
for typical Ref. d values (A = 3.0 a.u., i max = 6000 a.u.) 
resulted in a final error of 10%. 

An effective solution to the asymptotic growth prob- 
lem has been found, but first it is necessary to divulge 
further details of the trajectory implementation, as pre- 
sented schematically in Fig. [T] At t = 0, the two trajec- 
tory grids, i.e. "upper" (for ^ + ) and "lower" (for \&_) are 
coincident, with grid points at both xl and xr, and a uni- 
form grid spacing of Ax = {xr — xl)/(N— 1). Overtime, 
x]y, the right-most trajectory of the upper grid, evolves 
beyond the interaction region, as does xj - , the left-most 
point of the lower grid. At time £ s hift = Ax/w, when these 
extremal points have progressed one grid spacing beyond 
the interaction region, they are discarded, and new up- 
per/lower grid trajectories are introduced at xl/xr, and 
assigned ^± values in accord with Eq. |5]). The indexing 
for both grids is shifted by one, resulting once again in 
coincident grids identical to that at t = 0. This approach 
obviates the requirement of a nearest-neighbor search al- 
gorithm, thus greatly reducing CPU cost per time step. 
It does require that t s hift be a multiple of A, although in 
practice, this poses a very minor constraint. 

Note from Fig. [T] that at general times within the cy- 
cle described above, one or two points from each grid lie 
outside the range of the opposite grid. In such cases, 
one might consider resorting to extrapolation to evalu- 
ate the coupling term in Eq. although this approach 
is not recommended, 42 Instead, since the exterior points 
all lie near the edges of the interaction region, we sim- 
ply set V — > in Eq. ([3]). In principle, this results in 
unitary plane wave propagation (PWP) for the exterior 
points. In practice, the Euler time integrator used here is 
not unitary (nor is Runge-Kutta) — thus causing the ob- 
served growth of asymptotic probability. This situation 
is well understood in the literature, and various alterna- 
tive unitary time integrators exists However, since the 
analytical plane wave solution is known, we simply adopt 
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implementation is therefore as follows: 
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FIG. 1: Schematic indicating location of trajectory grid points 
(circles) over time, for upper (*!/+) and lower (4'—) grids of 
the constant velocity trajectory bipolar counterpropagating 
wave method. At t = 0, the TV = 5 points of both grids are 
distributed uniformly from xl to xr, with grid spacing Aa; = 
(xr — xl)/(N — 1). Over time, the upper/lower grid points 

move with constant velocity ±t) = ±y/2E/m. At time t s hiit = 
Ax/v, the two grids are again coincident in the interior region. 
The two outermost points are deleted, and replaced with two 
new points (indicated with open circles) at the opposite ends 
of their respective grids, to restore the situation at t = 0. 



the following exact unitary PWP scheme: 



*± fc (t + A) = tf ±fc (t)exp 



-EA 



(6) 



where 1 < k < N indexes individual trajectories. 

The above PWP modification, applied to extreme 
points on both sides of the two grids in lieu of Eulcr 
time integration, completely eradicates the asymptotic 
growth of probability over time. However, it leads to 
a new problem, i.e. asymptotic "ramping." The idea 
here is that as the x^ trajectory travels inward, calcula- 
tion of the plane wave portion of Eq. © shifts suddenly 
from the unitary PWP to the nonunitary Euler integra- 
tor, causing the final converged solution to manifest a 
gradual increase in asymptotic p+(x) from left to right. 
This spurious ramping was found to be on the order of 
a few percent for the test Eckart system — though the ef- 
fect is once again greatly exaggerated here, due to the 
first-order time integrator used. Still, even with realis- 
tic time integrators (e.g. fourth-order Runge Kutta) the 
ramping problem would cause difficulties for very high- 
accuracy calculations. A simple solution is nevertheless 
available: use the Eq. ^ PWP for the plane wave por- 
tion of Eq. ([3]) throughout the coordinate range (i.e. not 
just in the asymptotic extremes). The first-order forward 



V± k (t + A) =* ±fc (i)exp 



~EA 

Ti 



-lv(«fc)A[*±*(t) + * T (x^t)],(7) 

where \E , =F (x^,t) is computed via interpolation. Use of 
Eq. ([7} has eradicated the nonunitary error completely 
in the test Eckart application with parameters described 
above, resulting in a numerical solution accurate to sev- 
eral parts per thousand throughout the full x range. The 
error that remains adheres to the form of interpolation 
error, as discussed above. 



B. Interpolation 

As interpolation is a primary source of numerical error, 
it is essential for the present purpose that this operation 
be performed in as efficient and accurate a manner as 
possible. One of the key presumptions of successful in- 
terpolation is that the underlying mathematical function 
be smooth and well-behavedjSjiS, This results not only in 
more accurate computation, but also better characteri- 
zation of the numerical analysis. Fortunately, the bipo- 
lar CPWM methodology is designed to provide smooth, 
slowly- varying field functions even when *F itself is not; 
these have indeed been observed in all cases considered 
thus far, with the F constant-velocity functions being 
particularly slowly- varying. Other desired attributes of 
the interpolation routine include generality, robustness, 
systematic error reduction, and a minimum of required 
user intervention.— 

In the previous work^ fourth-order local polynomial 
interpolation was used, analogous to a moving least 
squares method with a stencil of five grid points for the 
interior of the grid.^2 At the grid edges, fewer sten- 
cil points (and lower-order polynomials) were used. In 
principle, if the field functions are sufficiently smooth, 
and the grid points ideally distributed, then higher-order 
calculations can lead to more accurate results. However, 
this is not guaranteed, and in some cases higher-order in- 
terpolation may introduce unwanted wild oscillation be- 
tween grid points.— As robustness and ease-of-use are key 
for the present purpose as described above, we simply 
adopt third-order interpolation as the most reasonable 
(and standard) choice. Also, since edge effects appear to 
be so important (Sec. Ill B| ). we retain third-order interpo- 
lation to the edges of the grids, i.e. four-point stencils are 
used throughout for polynomial interpolation. We also 
consider spline interpolation,-- for which the "stencil" 
is always two grid points, in effect, regardless of order. 
Spline interpolation has the advantage of continuous sec- 
ond derivatives throughout the coordinate range, yielding 
"stiffcr" interpolants that are more stable, and generally 
more accurate for smooth functions. 

The numerical implementation described in Sec. IIII Al 
employs interpolation only when evaluating the coupling 
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term, i.e. the very last \I/zp (xf , t) quantity in Eq. J?}- 
Although Eqs. © and ([7]) are written in terms of 
in fact it is the real- valued field quantities, r± and s±, 
obtained from 

®±(x) = r±(x)e ls±{x)/h , (8) 

that are designed to be most smooth and slowly-varying 
(especially in the classical limit).— Numerically, one 
might consider storing r±k and s±k values instead of 
^±fc's, but the corresponding Eq. ([7]) becomes compli- 
cated, nonlinear, and numerically less stable. 



Instead, we adopt a "best of both worlds" approach, 
whereby ^±k values are stored and used in Eqs. ^ and 
(O, r^'s and Szp/t's obtained at each time-step, and 
the latter interpolated to compute the coupling term in 
Eq. (O viaEq. ([8]). Obtaining r Tk = |* Tfc | is straightfor- 
ward; however, s^k = fa arctan (Im^I'i F fc/Re\E' = pfc) requires 
special care in the choice of branch cut, to ensure con- 
tinuity of the field function throughout the coordinate 
range. Our simple solution to this well-known dilemma 5 
is the recursive implementation 



Sip A: 

fa arctan 



Im4 



T(fc+1) 



Re* 

s^k + fa arctan 



Im(^ T(fc+ i)/W Tfc ) 



Re(*^ 



for * T (fc + i) = 
for * Tfe = 0; 

otherwise, 



(9) 



with ^zpo = by definition. Note also that due to the 
regular structure of the grid, there is no need for a search 
algorithm to compute the stencil for a given interpolant 
point, i.e. the total CPU effort per time-step associated 
with interpolation scales as TV rather than TV log TV, plac- 
ing it on a par with potential function evaluation. 

The above r/s interpolation scheme yields phenome- 
nal improvements in accuracy when applied to the test 
Eckart problem of Sec. IIII A[ as compared to the \l/± in- 
terpolation applied there. The numerical solutions still 
exhibit asymptotic oscillations; however, these are no 
longer constant in amplitude, but damp exponentially 
with increasing |a;|, as do the true mathematical solu- 
tions. In any case, at the — xl = Xr = 3.0 a.u. grid-edge 
values considered, the amplitude has been reduced from 
a few parts per thousand (Sec. IIII Ap to a few parts per 
million, using a grid of only N — 31 points. 

We conclude this subsection with a brief numerical 
comparison between polynomial and spline interpolation 
(both third order) in the context of the above algo- 
rithm as applied to the Eckart B problem with E = 
0.40Vb fSec. IIV A"|) . First though, as any complete spline 
implementation requires specification of the boundary 
conditions, we compared the accuracy of zero-second- 
derivative boundaries (natural splines) vs. the more 
usual zero-first-derivative case, for known functions sim- 
ilar to those computed below. The natural spline case 
was found to be about an order of magnitude more ac- 
curate throughout the coordinate range, and is therefore 
adopted. 

For comparative purposes, both natural spline and 
polynomial interpolation were applied to the Eckart B 
problem. Both cases are characterized by quick conver- 
gence with respect to grid size, TV, with the spline case 
converging more quickly. With respect to the time step, 



A, polynomial interpolation converges more quickly at 
first — i.e., for larger A values and lower accuracies. Be- 
yond a certain point however, the polynomial error flat- 
tens out and even increases with subsequent decrease in 
A (even if TV is allowed to increase), whereas the spline 
error becomes arbitrarily small. Further analysis reveals 
that this behavior of the polynomial interpolation is due 
to asympotic oscillations, which increase in amplitude 
with decreasing A, whereas in the spline case, the oscil- 
lation amplitudes are independent of A. The polynomial 
case thus exhibits an effective minimum achievable er- 
ror, on the order of 10~ 5 for the case considered (note: 
in conjunction with first-order forward time evolution). 
The parameters used for this optimized polynomial cal- 
culation were also applied to a spline calculation, which 
was found to be both more accurate (orders of magnitude 
so for the transmission probability) and faster (by about 
30%). Taking these facts into consideration, we adopt 
natural cubic spline interpolation from here on out. 



C. Discontinuous Constant Velocity Scheme 

The constant-velocity trajectory method as discussed 
above is not applicable for asymptotically asymmet- 
ric potentials 3 — i.e., Vl = ^ Vr, where Vl = 
linix—j-oo V(x) and Vr = lim^^oo V(x). Mathematically, 
this is because the coupling term in Eq. ([3]) does not 
vanish in the product asymptote (x — > oo), and so the 
ty±(x,t) components do not converge over time. Physi- 
cally, it is because the trajectory velocity v is not com- 
mensurate with the product translational kinetic energy, 
(E — Vr). In any event, the presumption Vl — Vr is 
too restrictive in practice to be of general use, and so 
some appropriate resolution must be found. One might 
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consider the classical trajectory bipolar CPWM, for in- 
stance, which would work fine for barrierless reactions; 
however, the more general and realistic case of an asym- 
metric reaction with a barrier (Fig. [3]) — and therefore 
below-barrier energies — would seem to be beyond the 
reach of either method. 

On the other hand, the above difficulty is simply a 
manifestation of a standard concern that is not at all 
new in reactive scattering theory — i.e., that the prod- 
uct and reactant asymptotic potentials, and associated 
quantum states, are not the same^ 3 - Various remedies 
have been prescribed for this, including use of a "divid- 
ing surface" 34 > 35 that partitions configuration space into 
reactant and product sides. For 1 DOF applications, the 
dividing "surface" is just a single point, labelled xq. We 
can easily incorporate the dividing point idea into the 
constant velocity scheme by using trajectories with reac- 
tant asymptotic velocities in the x < xq region, and an- 
other set with product asymptotic velocities in the x > xo 
region. The appropriate generalization of Eq. ([3]) (which 
holds for any asymptotic values of the potential) is 



L/R± 



dt 



jr{E-V- V L/R ) * L/R± 



i 

h 



v - V L/R ) *i 



AR=F> 



(10) 



with "L" referring to the reactant (left) side of xq, and 
"i?" the product (right) side. The corresponding trajec- 
tory velocities are v L / R — y/2(E — Vl/r)/™,. 

Although '5 itself, and its first spatial derivative, must 
be continuous across x = xq at all times, the components 
ty± are decidedly discontinuous. In a trajectory imple- 
mentation, ^l+^xq) and ^S> r ^(xq) are known quantities, 
and ^l-(xo) and ^ R+ (xo) are unknown; however, the 
^/ constraint above can be used to uniquely specify the 
latter quantities in terms of the former 4 For the sake of 
numerical efficiency, it is advisable that 1 $>l+ and ^r- 
trajectories cross xq at exactly the same time; when this 
occurs, these "edge" trajectories (for their respective re- 
gions) are destroyed, and new edge trajectories for 
and ^r+ simultaneously created without needing to re- 
sort to extrapolation. This is also when trajectory "shift- 
ing" occurs, as per Sec. lIII Al The above trajectory cross- 
ing condition requires that the grid spacings be different 
on either side of xq, i.e. Ax^/ R — ^ s hift^ L/_R- I n other 
respects, the algorithm is similar to Sec. IIII Al — except 
that when Xk ^ xq, we can not use Eq. ([6]) for the edge 
trajectories near xo, because the coupling in this region 
is significant. Instead, we must use Eq. ([7|), even though 
this implies extrapolation within a range (Axl/ r )/2 of 

Xq. 



D. Ramp Trajectory Scheme 

It would be beneficial to develop a single, robust, and 
efficient bipolar CPWM scheme with none of the limita- 
tions of the previously-described methods. The method 



should in other words be continuous, slowly-varying and 
extrapolation-free, yet applicable to completely general 
potentials such as that of Fig. [3j From the discussion 
in Sec. Ill Al it is natural to consider the generalized 
F approach in this regard. Specifically, the "semiclas- 
sical" trajectories considered need have no direct rela- 
tion to the actual potential V(x), but rather, are defined 
for some arbitrary effective potential, V e g(x), as v(x) = 
y/2 [E — V c s(x)} /m. The specific choice V g(x) — V(x) 
reproduces the classical trajectory scheme, V e s(x) = 
the constant-velocity scheme, and V e g(x) = Vl + {V R — 
Vi)Q{x — Xo) the discontinuous scheme. 
The ideal V e g(x) should: 

1. approach the flat V(x) function in both asymptotic 
limits, i.e. limx-^oo V e s(x) = V L / R 

2. vary smoothly and monotonically through the in- 
teraction region. 

Condition 1. results in vanishing asymptotic coupling, 
whereas 2. ensures there are no turning points, so that 
barrier tunneling poses no difficulties. A sigmoid ramp 
function is thus clearly indicated; below, we present a 
specific form based on the tanh function that is applica- 
ble to extremely general situations. Although this form 
might not be ideally suited to systems with reaction in- 
termediates, one could in that case add several ramp 
potentials together to match the reaction profile appro- 
priately. In any event, we must first work out the the- 
ory and numerical implementation for the generalized F- 
based CPWM approach. 

The generalized F decomposition of Eq. |lj for station- 
ary scattering states in 1 DOF is uniquely determined 
from the time-independent Schrodingcr equation and the 
relation 3 -^ 



Zv n 



(11) 



For completely general v(x), these lead to 



where 



I 1 I 1 A 

± —mv ± —A 

h 2 



4-±±-A* T , (12) 



v cS ) 



(13) 

For consistency, it is convenient to reexpress the v 
derivatives above in terms of V e g and its derivatives. If 
the trajectory equations of motion are to have no ^± 
spatial derivatives as desired, then we must construct a 
convective term somehow; this is easily accomplished by 
multiplying Eq. (fT2"| by ±i>. The result is: 



± 1 -v 
4 \E-V eS 



+ % - (mv 2 - V + V eS + C) 



[y-T4ff-C]* T , where (14) 



8 



2m 



■5 

Tg 



E-V r , 



V" 



E-V R < 



in Sec. Ill At we take 
v / solution, i.e. 



0) to be the basic WKB 



Next, we replace mv 2 in Eq. (fT4|) with 2(S— y c ff). Adding 
— {i/%)Ety± =F t>'J , ' ± to both sides, and recognizing that 
for the stationary solution, d^>±/dt = — (i/h)E^/±, we 
obtain finally 



~dt 



4(^fc) + ^ ( ^- y - yeff+a) 



Note that for constant velocity trajectories, C = V c s = 
0, and Eq. (TT6"|) reduces to Eq. (g]), as it should. Also in- 
teresting (and useful) is the fact that Eq. (TO))) satisfies 
the combined continuity relation— — i.e. the total com- 
bined probability (p+ + P-) is conserved over time, even 
though the individual p± are not. The precise flux rela- 
tion giving rise to this is 

^ = -J± ± ~ (V - V e g - C) Im [* + * , (17) 

where j± — ±vp± is the flux, defined in terms of the gen- 
eralized trajectory velocities. For the constant velocity 
case, this reduces to Ref. 0Eq. (13). 

In implementing the above scheme numerically, we bor- 
row from the successful constant-velocity algorithm as 
much as possible. The condition that over a duration of 
time, tshiftj each grid point should progress to the pre- 
cise spot vacated by its nearest neighbor, implies a very 
nonuniform grid — with a higher density of points in areas 
where v is lower, as is appropriate. Though by no means 
necessary, we find it convenient to compute all of the 
trajectories in advance. Actually, only the single trajec- 
tory x± (t) need be computed — starting at x = Xz, and 
propagating classically in accord with V e g to a bit past 
x — xji. The other trajectories can then be obtained 
via x^(t) = x~\ [t + (k - l)t s hift], and x^(t) = x^(—t). 
The xf (t) values for each positive integer multiple of the 
CPWM time step, A, are computed and stored during 
this preprocessing phase, as are the corresponding effec- 
tive potential values V e g, V^ ff , and — a technique that 
can in principle save much CPU time during the subse- 
quent CPWM propagation phase. Note that by design, 
the trajectory version of Eq. (|16p involves no explicit spa- 
tial differentiation of fy±; once again, the largest source of 
numerical error (and CPU effort) is interpolation. To en- 
sure that the trajectory calculation does not contribute 
appreciably to the error, it is performed using a much 
smaller time step, <5 = A/100; even so, the CPU cost of 
the preprocessing phase is a trivial part of the whole cost. 
Also, numerical comparisons with S = A/30 and A/300 
reveal no discernible difference in computed results. 

Another difference from the constant-velocity scheme 
is the choice of initial conditions. Based on the discussion 



* + (x, 0) = y/v L /v{x)exp [is (x)/h] , (18) 

where the initial action, sq(x), is also computed during 
the preprocessing phase (the initial reflected wave is of 
course still zero). The generalization of Eq. ([7]) is seem- 
ingly non-trivial: 



* ±fc (i + A) = * ±fc (t) 



(16) 



(19) 



exH-A 



E- 2U Gff ± -hv 



■-A(V- Veff" 

n 



with all functions evaluated at 



4"' \E-V eS 
C) [* ±k (t) + ^(xf,t)] 

f. The rationale is the 
same as before however: the asymptotically vanishing 
coupling term in Eq. (fT6|) (but for 'J', rather than just 
is treated using forward Euler propagation; everything 
else is exponentiated as per Eq. (|6|). 

The proceeding discussion is completely general, i.e. 
applies to any choice of V e g(x). However, for single- 
barrier reactions (no intermediates), we advocate the fol- 
lowing analytical ramp functional form, whose deriva- 
tives are also specified: 

V e g(x) = ( J R ~ Vl ) {tanh [(3(x - x )} + 1} (20) 



( Vr-Vl 
V 2 J {cosh 2 [f3(x-x )} 



I cosh [p(x — xo)\ 



As in Sec. IIII C[ xo should be in the middle of the V(x) 
interaction region, and 1//3 comparable to the interaction 
region size. 

E. Summary and Numerical Recipes 

To summarize, the key numerical developments in this 
section have been: (1) use of trajectory grids that are 
coincident when t is a multiple £ s hift (implying that £ s hift 
itself is a multiple of A); (2) unitary plane wave (or ex- 
ponential) propagation throughout the coordinate range; 
(3) third-order spline interpolation of r± and s±, to 
evaluate the coupling contribution to the time-evolution 
equations. 

For asymptotically symmetric potentials, the recipe for 
numerical propagation of the constant velocity trajectory 
quantities is as follows: 

• Move each grid point, x k , via xf (t + A) = x J (t) ± 
vA 

• For exterior grid points outside the range of the 
opposite grid, update ^>±k via Eq. ©. 
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• For interior grid points, update ^>±k via Eqs. 0, 
©, and (9). 

For asymptotically asymmetric potentials, the ramp 
trajectory method is preferred over the discontinuous tra- 
jectory method (Sec. lIV( l. The ramp trajectories are com- 
puted and stored a priori, in a preprocessing step. The 
subsequent numerical propagation is as follows: 

• For exterior grid points outside the range of the op- 
posite grid, update ^±k via the first line of Eq. (|2"0)) . 

• For interior grid points, update ^± k via Eqs. (f^Uj) . 
©, and (9). 



IV. RESULTS 

In this section, the three different numerical algorithms 
described in detail in Sec. IIIH — i.e. the constant veloc- 
ity, discontinuous constant velocity, and ramp trajectory 
versions of the non-wave-front F trajectory method — are 
applied to a variety of test applications, in order to assess 
them for robustness and numerical efficiency. In all cases, 
first-order forward Euler time-evolution is employed as 
discussed previously, leading to "artificially" small time 
steps, A (e.g. as compared to the fourth-order Runge- 
Kutta results of Ref. y). The mass m = 2000 a.u. is used 
in all applications. 

A thorough and detailed convergence study was per- 
formed for each application as follows. First, given the 
small oscillations with x that characterize the numeri- 
cal p±(x) solutions in the asymptotic regions, the usual 
determination of reflection and transmission probabili- 
ties as P ro fi = |*_(X L )| 2 ; Ptrans = {v R / V L ) + {x R )\ 2 

is not used. Instead, a small range of points near xl/r 
is consulted to characterize both mean values and un- 
certainties for these quantities. Reducing the oscillatory 
uncertainty to an acceptably small level is one of two con- 
ditions that must be satisfied to achieve convergence; the 
second is that variation of the mean P rcr i and Ptrans values 
with respect to all numerical parameters be acceptably 
small. The convergence parameters, in the cyclical order 
in which they are varied (keeping all others fixed) are 
as follows: maximum time, i max ; grid spacing, Ax; time 
step, A; grid boundaries, Xl/r- For each of the appli- 
cations considered, tens of convergence calculations were 
performed, usually over multiple cycles. 



A. The Eckart Barrier 

The first application considered is the symmetric 
Eckart barrier ) 44 i 45 

V(x) = Vbsech(ax) 2 , (21) 

with parameter values Vq = 400 cm -1 , and a = 3.0 a.u, 
which for the present paper is labeled "Eckart A." This 
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FIG. 2: Potential energy curves for several symmetric poten- 
tial systems considered in this paper: Eckart A barrier (solid); 
Eckart B barrier (dashed); double barrier (dotted). All units 
are atomic units. 

application is identical to the Eckart problem of Ref. d, 
and is chosen here in part to compare with the previous 
results. A plot of this potential is presented in Fig. [2] 

A detailed convergence study was performed for the en- 
ergy equal to the barrier peak, i.e. E = Vq = 400 cm -1 s» 

0. 0018 hartree. This energy was chosen because it pro- 
vides a good splitting of probability between P rcn and 
Ptrans- Of the three numerical methods considered, only 
the constant velocity method is applicable here — because 
both the discontinuous constant velocity and ramp tra- 
jectory methods reduce to the constant velocity method 
when the potential V(x) is symmetric. In any event, the 
parameters that resulted from the convergence procedure 
described above, with 10 -4 as the target accuracy for 
both Pcfi and Ptrans, are presented in TableJI] column 3. 
These can be compared directly to the Ref. calculation 
(computed to the same level of accuracy), presented in 
column 2, as well as to the analytical values given in the 
last two rows. 

The table indicates that the new algorithm requires 
substantially fewer grid points than the old — a testa- 
ment to the efficacy of the new PWP and interpolation 
routines. However, the new algorithm requires a much 
smaller time step (more than 50x), owing to the fact that 
the non-PWP time integrator is only first-order, rather 
than fourth-order. Presumably for the same reason, A 
was found to be the slowest parameter to converge. From 
Table HI it is also clear that the computed P re fi value is 
substantially more accurate than Ptrans- Since only one of 
these two quantities is needed in practice (because their 
sum is unity), we performed a second convergence study 
for which P rcn alone was converg ed to 10" 4 . The results, 
presented in Table U column 4, indicate a very substan- 
tial efficiency improvement, requiring a mere N — 13 grid 
points and 0.88 seconds on a 2.60 GHz Pentium CPU — as 
compared to N = 31 and 5.0 seconds for column 2. Note 
that for both calculations, the oscillatory error is sub- 
stantially smaller than the overall convergence error — 

1. e. it would have been essentially just as accurate to use 
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TABLE I: Convergence parameters for constant velocity trajectory method applied to several symmetric potential systems, 
with energies near or at the barrier height. All units are atomic units. Column 2: older algorithm used in Ref. 3. Columns 3-6: 
newer algorithm of Sec. IIII Al Convergence of quantities indicated is to 10 -4 for Eckart A calculations, and 10 -3 for others. 
Rows 7 and 8: computed reflection and transmission probabilities, including oscillatory uncertainties. Last two rows: exact 
reflection and transmission probabilities for Eckart systems. 



On Antitv 






Symmetric potential system 




and 




Eckart A 




Eckart B 


double Gciussiciii 


symbol 


Ref. 3, 


Prcfl and Ptrans 


Prefl only 


P rc fl only 


Prefl and Ptrans 


grid size, N 


31 


20 


13 


25 


20 


left edge, xl 


-3.0 


-2.0 


-2.0 


-3.0 


-3.0 


right edge, xr 


3.0 


2.0 


1.5 


2.1 


2.5 


grid spacing, Ax 


0.20 


0.211 


0.292 


0.213 


0.289 


time step, A 


10.0 


0.156 


0.324 


0.046 


0.612 


max time, t max 


10000 


3899 


4105 


3204 


41015 


computed P re fl 




.28336 ± .00004 


.28323 ± .00010 


.4587 ± .0002 


.7936 ± .0002 


Computed Ptrans 




.71646 ± .00005 


.71625 ± .00005 


.5355 ± .0003 


.20498 ± .00001 


exact P re fl 




0.283358 


0.283358 


0.459605 




exact Ptrans 




0.716642 


0.716642 


0.540395 





-Ptrans = P+{xr), etc. 

Finally, column 5 of Table [TJ presents results for the 
much broader and taller "Eckart B" barrier, i.e. with 
parameter values Vq = 0.011 hartree (« 2414cm -1 ), and 
a = 1.364 a.u. The Eckart B potential is also presented 
in Fig. [21 The goal here is to determine whether the 
new algorithm is robust across a wide range of system 
parameters, and to what extent the much larger-action 
Eckart B system requires additional computational ef- 
fort. The Eckart B system is also used to analyze deep 
tunneling performance in Sec. IIVB1 Here, though, we 
consider only the barrier energy E — Vq. A comparison 
between columns 4 and 5 reveals that a larger grid and 
smaller time step are needed in the Eckart B 
would be expected based on system parameters. How- 
ever, the near order-of-magnitude time step reduction 
is far greater than expected — particular given that the 
Eckart B calculation is less accurate to one digit of pre- 
cision. This is due to the convergence of computed error 
with decreasing A, which, though initially fast, is found 
to become quite slow (essentially inverse scaling) beyond 
a certain threshold accuracy level. For Eckart B, this 
threshold appears to be around 3 x 10~ 3 . In any case, 
higher-order time evolution is expected to resolve this 
difficulty completely. 



B. The Deep Tunneling Regime 

To analyze the efficacy of the constant velocity method 
in the deep tunneling regime, the Eckart B problem was 
solved for a range of below-barrier energies, E < Vq. The 
large size of this barrier ensures that even energies rela- 
tively close to the top of the barrier — e.g., E = 0.8 Vq — 
are characterized by quite small transmission probabil- 
ities, Ptrans- By E = 0.1 Vq, the Ptrans values are ex- 
tremely small (on the order of 10 -9 ). Detailed conver- 
gence studies were performed for the above two energy 



values, as well as for the intermediate value E — 0.4 Vo, 
in order to assess the method across the whole energetic 
range. The details are presented in Table [Til For sim- 
plicity, xl = xr was presumed. The quantity Ptrans was 
converged to an absolute accuracy that varied substan- 
tially with E, as indicated in Table ITTI row 6. In relative 
accuracy terms, this translates to three-to-five significant 
digits for each computed Ptrans value, with the E = 0.4 Vo 
case substantially more accurate than the other two. In 
all cases, the convergence error is comparable to the ac- 
tual error, obtained relative to the exact Ptrans values 
presented in row 8. These findings are particularly satis- 
fying for the extreme E = 0.1 Vq case, which is computed 
here to an absolute error of 10~ 12 or so. 

The convergence trends with respect to decreasing E, 
as evident in Table (TTJ are illuminating. Perhaps most 
striking is the fact that only a modest increase in co- 
ordinate range, xl/r, is observed — completely unlike 
the usual situation for conventional quantum scattering 
methods (see below). A modest increase in grid size, 
N, arising primarily from a reduction in grid spacing, 
Air, is nevertheless required. This is true, even though 
the wavelength increases with decreasing E, because the 
lower energy calculations must be performed to a sub- 
stantially higher level of absolute accuracy. Comparing 
rows 1 and 6, the convergence appears to be exponentially 
fast, with something like twelve extra points needed for 
each additional digit of accuracy. Relative accuracy also 
plays a role in this regard, however, e.g. the 0.8 Vq calcu- 
lation done to the same absolute accuracy as the 0.4 Vo 
calculation is found (not presented in Table [TTJ) to require 
a greater number of grid points than the latter. 

There is a very substantial growth of i max with decreas- 
ing energy, owing to the fact that the trajectories move 
much more slowly. However, for the same reason, there 
is also a substantial increase in the time step size, A, so 
that the growth in the total number of time steps (and 
CPU effort) is much more modest. The above comments 
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TABLE II: Convergence parameters for constant velocity trajectory method applied to Eckart B system over a wide range of 
below-barrier energies, E < Vq. All units are atomic units. Convergence is of transmission probability, P tr ans, only, to accuracy 
level indicated in row 6. Computed Ptrans values including oscillatory uncertainty are listed in row 7, with exact values given 
in row 8. 



Quantity 
and symbol 




Energy E as fraction of barrier height Vq. 




0.8 U 


0.4 Vq 


0.1 Vq 


grid size, N 


33 


88 


133 


grid edges, x L/R 


T3.0 


T4.0 


T4.0 


grid spacing, A:r 


0.188 


0.092 


0.061 


time step, A 


0.126 


0.097 


1.156 


max time, f ma x 


4424 


46424 


462285 


target accuracy for Ptrans 


0.02(-2) 


0.0005(-5) 


0.01(-10) 


Computed Ptrans 


4.444±.020(-2) 


1.5588±.0004(-5) 


9.920±.002(-10) 


exact Ptrans 


4.462(-2) 


1.5594(-5) 


9.920(-10) 



apply for calculations at comparable relative accuracies, 
e.g. Table HJ columns 2 and 4. The intermediate calcu- 
lation (column 3), done to a substantially higher relative 
accuracy than the others, is found to be in the regime 
of slow convergence with respect to A (Sec. IIV Ap and 
therefore manifests a much smaller converged value for 
this parameter. 

The extreme E = 0.1 Vq calculation represents a very 
challenging test case for conventional quantum dynamics 
methods. For comparative purposes, an optimized sinc- 
DVR calculatio n 34 ! 35 ' 46 with quartic absorbing bound- 
ary conditions was also performed. Such methods al- 
ways increase the coordinate range beyond the interac- 
tion region — to a degree proportional to the wavelength, 
and thus very substantial at low energies. For instance, 
x l/r — T40.0 a.u. for the DVR calculation, a full or- 
der of magnitude larger than for the CPWM calculation. 
The resultant grid size, N — 700, was also much larger — 
particularly given that CPU effort scales as N 3 for DVR 
methods, but only as N for the present approach. 

The above discrepancies are magnified even more for 
shorter, narrower barriers than Eckart B. In this context, 
deep tunneling is observed only when E is a tiny fraction 
of V a . For the Eckart A problem with E = 10~ 4 U , 
for instance, attempts to perform a converged sinc-DVR 
calculation proved unsuccessful, owing to the enormous 
coordinate range required. Instead, an approximate, 
semiclassical tunneling calculation was performed, giv- 
ing rise to the prediction Ptrans = 3.82 x 10 -3 . In fact, 
this is more than lOOx larger than the actual value, 
-Ptrans = 2.851 x 10~ 5 . These examples serve to under- 
score the difficulty of computing accurate deep tunnel- 
ing transmission probabilities using conventional meth- 
ods, which may be greatly alleviated with the present 
bipolar CPWM approach. 



C. The Uphill Ramp 

The first asympotically asymmetric system consid- 
ered here is the continuous "uphill ramp,"— defined 
via Eq. ([2U)l and the parameter values Vl = 0, Vr = 
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0.00J - 
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FIG. 3: Potential energy curves for two asymmetric poten- 
tial systems considered in this paper: uphill ramp (solid); 
barrier ramp (dashed). Dotted line indicates energy value, 
E = 0.0023, used for both systems. All units are atomic 
units. For the ramp trajectory calculations for both systems, 
the trajectory potential, V c g(x), is the solid curve. 

400 cm- 1 w 0.0018, x = 0, and [3 = 2.5 a.u. Thus, 
V(x) — V e g(x), so that the ramp trajectories are them- 
selves classical trajectories. The energy E = 0.0023 a.u 
« 500 cm- 1 is deliberately chosen to be only moderately 
larger than Vr 7 to accentuate the asymmetric difference 
between Vl and Vr. A plot is presented in Fig. [3] The 
system and energy parameters are identical to those con- 
sidered in Ref. y, and are also comparable to those of 
the Eckart A calculation in Sec. IIV Al enabling a fairly 
direct comparison. Detailed convergence studies were 
performed for both the discontinuous constant velocity 
and ramp trajectory methods, as presented in Table \\T\\ 
Both P ren and Ptrans were converged, with the former to 
a higher level of accuracy, appropriate to the different 
scales of the two quantities. Note that this is the first 
example considered for which Ptrans 3> -Prcfl- 

The discontinuous constant velocity case is considered 
first, i.e. Table IIIII column 2. As compared with the 
Eckart A calculation in Table Q] column 3, the most com- 
pelling difference is that a much higher density of grid 
points is needed (N = 80 vs. N = 20), even though the 
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TABLE III: Convergence parameters for discontinuous constant velocity and ramp trajectory methods, applied to two asym- 
metric potential systems with energy E = 0.0023 a.u ~ 500 cm -1 slightly larger than asymptotic potential difference 
(Vr — Vl) = 400 cm -1 . All units are atomic units. Columns 2 and 3: uphill ramp system. Columns 4 and 5 barrier 
ramp system. Convergence is of both P rer i and Ptrans to accuracy levels indicated in rows 7 and 8. Last two rows: computed 
reflection and transmission probabilities, including oscillatory uncertainties. 



Quantity 




Asymmetric potential system 




and 


uphill ramp 


barrier 


ramp 


winhnl 


discontinuous 


ramp traj. 


discontinuous 


rftmn tra i 

i din yj ur < ' I ■ 


grid size, N 


80 


25 


80 


23 


left edge, xl 


-2.0 


-2.0 


-2.0 


-2.0 


right edge, xr 


2.0 


2.0 


2.0 


2.5 


left grid spacing, Axl 


0.080 


0.252 


0.080 


0.315 


time step, A 


0.235 


0.500 


0.754 


0.300 


max time, £ max 


6594 


5478 


7912 


7474 


target accuracy for P re fl 


.00004 


.00002 


.0002 


.00008 


target accuracy for Ptrans 


.0004 


.00015 


.0006 


.0004 


computed P ren 


.023838 ± .000002 


.023919 ± .000003 


.45452 ± .00003 


.45454 ± .00003 


Computed Ptrans 


.97557 ± .00001 


.9762 ± .0001 


.54465 ± .00003 


.54602 ± .0001 



final computed accuracy is somewhat lower. This is most 
likely attributable to the use of extrapolation (Sec. IIII C|) 
in the middle of the interaction region where the coupling 
is greatest. One indication is the fact that the asymptotic 
oscillations are much smaller than for Eckart A, and not 
a contributing factor to the total error. Another signifi- 
cant difference is found in the parameter i max , which is 
larger in the uphill ramp case, probably owing to the fact 
that product trajectories are substantially slower than for 
Eckart A. Perhaps for the same reason, A is also some- 
what larger for the uphill ramp. 

The ramp trajectory convergence parameters are 
shown in Table IIIII column 3. The N = 25 grid size is 
now greatly reduced in comparison to column 2, almost 
to the Eckart A level, again suggesting that the primary 
source of error in the discontinuous case is extrapola- 
tion. The small increase in N relative to Eckart A may 
be due to the fact that grid spacing in product and reac- 
tant regions are not independently adjustable (which may 
also be a contributing factor in the discontinuous calcula- 
tion). Another key difference between the two asymmet- 
ric calculations is that i max is substantially reduced in 
the ramp trajectory case — most likely reflecting the fact 
that the initial semiclassical approximation \l/+(x, 0) is 
much closer to the actual solution (Sec. IIIID[) . Never- 
theless, the time step A is substantially larger. Taken 
together, all of the above findings imply that CPU effort 
for the ramp trajectory calculation is greatly reduced in 
comparison to the discontinous calculation. 



D. The Barrier Ramp 

The barrier ramp system (Fig. [3]) is the first asymmet- 
ric reaction with a barrier to be treated using CPWM 
methods. Serving as representative of a generic reaction 
profile, moreover, it is an extremely important bench- 
mark system. This is especially true for the ramp trajec- 



tory method, because V c s(x), being monotonic, can no 
longer be equal to V(x), which might in principle affect 
computational efficiency. To construct a suitable barrier 
ramp potential, we simply add an Eckart barrier and an 
uphill ramp together. The Eckart parameters [Eq. (|2ip ] 
are taken to be Vq = 0.0015 a.u. and a — 2.5 a.u, 
whereas the uphill ramp is identical to that of Sec. IIV CI 
The ramp trajectories are also identical to Sec. IIV CI i.e. 
V e s(x) is the same as before. The energy E = 0.0023 a.u 
is also the same as in Sec. IIV CI as are the asymptotic 
trajectory velocities, w_L/fl- 

The discontinuous and ramp trajectory convergence 
parameters are shown in Table IIIII columns 4 and 5, 
respectively. In both cases the grid parameters are 
nearly identical to the corresponding uphill ramp calcu- 
lations (columns 2 and 3) — implying that the presence of 
the barrier introduces no numerical difficulties for cither 
method. The i ma x values are somewhat larger, undoubt- 
edly owing to the fact that the barrier ramp case has 
much greater reflection. 2 ' 3 For the discontinuous calcu- 
lation, A is lower than the uphill ramp value, perhaps 
because the absolute accuracy of the computed P rcu is 
much less. This situation is reversed for the ramp tra- 
jectory calculation, for which a somewhat smaller A is 
required. The reason may have to do with the fact that 
V e s(x) ^ V(x), giving rise to larger coupling in Eq. (fl"6j) . 



E. The Double Barrier 

The last case to be considered is a system with reac- 
tion intermediates. In Ref. [3], a double-Gaussian barrier 
was employed. Here however, we wish to have more flex- 
ibility in the functional form, in order that a wider range 
of reaction profiles may be considered. To achieve this, 
a sum of two barrier ramps is used, centered at the tran- 
sition states on either side of the reaction intermediate. 
Though symmetry is not required, for simplicity we adopt 
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the symmetric form here as follows: 
V(x) = V {sech[a(x + x )] 2 + sech[a(x ~- x a )} 2 } + 
(Va/2) {tanh \j3(x + x )] - tanh \fi(x - x )}} 

(22) 

The parameters Vo = 0.0015 a.u, x — 1.0, a — (3 = 2.5 
a.u., and Va = 0.0005 a.u. are comparable to those used 
for the barrier ramp and Eckart A systems. The energy, 
E = 0.0014 a.u, lies slightly below the barrier peak of 
around 0.00175 a.u. A plot is presented in Fig.[5J 

Although this system would be a good candidate for 
a double ramp trajectory scheme (Sec. MI D[) . the con- 
stant velocity trajectory algorithm is used here. Con- 
vergence for both P re fi and -Ptrans to 10~ 3 is indicated 
in Table U column 6. This accuracy is somewhat lower 
than for Eckart A (column 3), to avoid slow A conver- 
gence (Sec. IIV Ap . The coordinate range is larger than 
for Eckart A, but the grid density is actually somewhat 
smaller (despite the fact that the potential has more fea- 
tures), such that the N — 20 grid size is the same for 
both. The convergence time £ max is much larger than for 
Eckart A, as is expected for a system with tunneling and 
reaction intermediates^ In other respects, the calcula- 
tion is similar to, and no more difficult than, the other 
calculations without reaction intermediates. 

V. CONCLUSION 

Based on the results as presented above, we conclude 
that all of the primary goals of this paper have been 
achieved. In particular, the robustness of the methods 
are demonstrated by the very broad range of test cases 
considered. These may be definitively regarded as "rep- 
resentative" of many other 1 DOF calculations (not de- 
scribed here) that were also performed. The efficiency 
of the methods is also established, in that most of the 
calculations that do not involve deep tunneling or reac- 
tion intermediates require no more than about one sec- 
ond of CPU time. Finally, the accuracy levels that can 
be achieved have been greatly expanded beyond previ- 
ous limits, e.g., to twelve-to-thirteen digits past the dec- 
imal in the extreme deep tunneling case (Table |TT] col- 
umn 4) — probably the highest accuracy ever achieved in a 
trajectory-based quantum calculation. Although 1 DOF 
systems are usually regarded as trivially soluble, these 
can nevertheless pose severe numerical difficulties when 
extremely deep tunneling is involved; it is therefore par- 
ticularly encouraging to see that the new methods work 
well in both this limit, and for more typical molecular 
situations. 

Indeed, there do not appear to be any 1 DOF appli- 
cations for which the present algorithms would not work 
well. This is of the utmost importance for the main mo- 
tivation of the present work, i.e. creating a parameter- 
free (except for convergence parameters) "black-box" 
quantum dynamics code for solving any 1 DOF station- 
ary scattering problem. A single universal algorithm 



is sought, although three specific implementations were 
considered in this paper: (1) constant velocity ; (2) dis- 
continuous constant velocity; (3) ramp trajectory. As 
both (2) and (3) reduce to (1) in the asymptotically sym- 
metric potential case, the only real choice is between (2) 
and (3). For reasons elucidated in Sec. IIV1 (3) is clearly 
the better choice, both with respect to numerical effi- 
ciency and robustness. As evidence of the latter, about 
the same number of grid points is required for systems 
with similar parameters, regardless of whether the poten- 
tial is asymmetric, or has one or more barriers (compare 
Table U columns 3 and 6 with Table IIIII columns 3 and 
5). 

The ramp trajectory scheme, (3), is therefore proposed 
for the universal purpose described above. One can imag- 
ine a straightforward procedure to automatize the selec- 
tion of suitable ramp parameters for any given potential 
without reaction intermediates. Such a procedure would 
presumably be easily generalized to multiple ramp tra- 
jectories in the case of reaction intermediates, thus main- 
taining "black-boxness" for all cases. On the other hand, 
even a single ramp is seen to perform very well for mul- 
tiple barrier systems [as in Sec. IIV El due perhaps to the 
smoothness of 14ff(x)], so that a single ramp may be all 
that is ever actually needed in practice. In fact an earlier 
preliminary investigation (without the numerical refine- 
ments of Sec. IIII[) obtained slightly faster convergence 
and smoother field functions for single vs. double ramp 
calculations. Note that the dynamical equations derived 
in Sec. lIIIDl are completely general, and could in princi- 
ple be used to explore other candidate functional forms 
for V eS (x). 

In any case, well-documented stand-alone fortran 
codes for all three algorithms considered here are avail- 
able from the author on demand — which, in the case of 
(3) at least, will also incorporate future developments. 
Note that the weakest algorithmic aspect at present is the 
first-order forward time integration of the coupling term 
in the dynamical equations — leading to slow A conver- 
gence beyond a certain point, and much smaller typical 
A values than in Ref. y (10 a.u.), which used standard 
fourth-order Runge-Kutta integration. The first order 
of business is therefore to incorporate a similar time- 
integrator into algorithm (3) — or perhaps something even 
more advanced, such as the time-adaptive Cash-Carp 
method^ or multistep Adams-Bashforth method^ Pre- 
liminary investigations indicate that this will lead to far 
greater time step sizes even than those used in Ref. H — 
e.g., A = 100 a.u. or more. If not the case already, such 
a modification would result in the most accurate and ef- 
ficient quantum dynamics codes available for 1 DOF sys- 
tems. 

Of course, the primary challenge for future develop- 
ment is multidimensional systems. As discussed in Sec.Hl 
the theoretical developments are already in place, both 
for stationary state and wavepacket scattering dynam- 
ics. These will be presented in detail in future publica- 
tions, but a brief outline may be found in the Summary 
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and Conclusions of Ref. 0. For the multidimensional sta- 
tionary state method, numerical instability of the sort 
described here has proven to be an issue. Accordingly, 
new multidimensional algorithms will be developed that 
incorporate the modifications of the present work. In 
addition to (hopefully) providing numerical stability, the 
reduced number of grid points der DOF that character- 
izes the present approach should lead to far greater re- 
ductions in CPU effort for applications with more DOFs. 
For wavepacket dynamics, it turns out that more conven- 
tional Bohmian trajectories are appropriate (defined via 
mv± — s'_j_). In anticipation of this, we have analyzed 
s' ± curves for the stationary solutions obtained here, and 
found: (1) smooth, but not always monotonic, behavior; 
(2) asymptotic agreement with trajectories used here for 



all but the ^_(x — > oo) asympotes; (3) exactly oppo- 
site trajectories in the ^-(rr — > oo) asymptote. The last 
observation, though curious, may not be significant, as 
r-(x) also vanishes in the same limit. 
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